(* ::Package:: *)

SetDirectory["C:\\Projects\\branches\\HMC_core\\output"]


param=OpenRead["..\\input\\params_hmc.txt"];
Skip[param,{Number}];
Tstart=Read[param,Number];
Tstep=Read[param,Number];
Tend=Read[param,Number];
Pacc=Read[param,Number];
Nmd=Read[param,Number];
Nmc=Read[param,Number];
Skip[param,{Number}];
Nruns=Read[param,Number];
Skip[param,{Number}];

Close[param];

temperatures=IntegerPart[(Tend-Tstart)/(Tstep)]+1;

templist=Range[Tstart,Tend,Tstep];

funcname="Griewank";
NmcToEquil=100;

sdOfEpotEnd=Array[c,{temperatures,Nruns}];
sdFinal=Array[b,{temperatures,2}];


If[Nmc<=NmcToEquil,Print["Nmc<=NmcToEquil in code"];Abort[];,Null];


index=1;
For[temp=Tstart,temp<=Tend,temp+=Tstep,
 Print["Index = ",index, "\nTemp = ", temp];
 For[i=1,i<=Nruns,i++,
  loadstring=".\\"<>funcname<>"\\HMC_diagnosis"<>ToString[Nruns*(index-1)+i]<>".out";

  fid=OpenRead[loadstring];

  Skip[fid,String];

  xsum=0.0;
  xsqsum=0.0;
  For[iNmc=1,iNmc<=Nmc,iNmc++,
   Skip[fid,{Number,Number}];
   If[iNmc>NmcToEquil,
    x=Read[fid,Number];
    xsum+=x;
    xsqsum+=x^2;,
    Skip[fid,Number];
   ];

   Skip[fid,{Number,Number,Number,Number}];
  ];

  Close[fid];

  Num=Nmc-NmcToEquil;

  sdOfEpotEnd[[index,i]]=Sqrt[(xsqsum/Num)-((xsum/Num)^2)];

 ];

 sdFinal[[index,2]]=Sum[sdOfEpotEnd[[index,j]],{j,1,Nruns}]/Nruns;
 sdFinal[[index,1]]=templist[[index]];
 index+=1;
]


ListPlot[Tooltip[sdFinal],PlotRange->{{Tstart,Tend},Full},AxesLabel->{Temperature,SD},AxesOrigin->{Tstart,0},PlotStyle->PointSize[Medium]]
